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ABSTRACT 

We describe and demonstrate a method for increasing the resolution locally in a 
Smoothed Particle Hydrodynamic (SPH) simulation, by splitting particles. We show 
that in simulations of self-gravitating collapse (of the sort which arc presumed to occur 
in star formation) the method is stable, and affords great savings in computer time 
and memory. When applied to the standard Boss & Bodenheimer test - which has 
been shown to depend critically on fulfilment of the Jeans Condition - the results are 
comparable both with those obtained using Adaptive Mesh Refinement, and with those 
obtained using a standard high-resolution SPH simulation, but they are achieved with 
considerably less computational resource. Further development and testing is required 
before the method can safely be applied to more general flows. 
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1 INTRODUCTION 

Numerical simulations of star formation — even the highly 
idealized simulations to which contemporary astrophysics is 
limited — require very large computational resources. The 
dynamics is self-gravitating, so the gravitational field has 
to be recalculated at each time-step. The ranges of density 



(~ 10 



g cm ^ to ' 



10^ g cm ^) and linear scale (~ 10^ 



cm to ~ 10 cm) are very large, so high resolution is needed. 
The geometry is complex, so three-dimensional simulations 
are essential. The dynamics is affected by a variety of radia- 
tive, thermal, chemical and magnetic effects, with the result 
that the energy equation is not a local function of state. 
The initial and boundary conditions are chaotic and poorly 
constrained by observation. 

Smoothed Particle Hydrodynamics (SPH) is a particle- 
based hydrodynamic scheme which accommodates some of 
these requirements automatically, by virtue of being La- 
grangian, having no imposed geometrical constraints, and 
being readily combined with particle-based gravity solvers, 
for example Tree Code Gravity (TCG; Barnes & Hut 1986; 
Hernquist 1987; Hernquist & Katz 1989). SPH is also rel- 
atively straightforward to implement, at least in its most 
basic formulation. More refined versions have to be used to 
obtain realistic results when there are large gradients of den- 
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sity (e.g. shocks and ionization fronts) and/or velocity (e.g. 
shocks and rapidly shearing discs). 

Since its initial realization (Lucy 1977; Gingold & Mon- 
aghan 1977), SPH has been developed and refined by several 
workers, for example Lattanzio et al. (1985; artificial viscos- 
ity). Nelson & Papaloizou (1993, 1994; variable h), Bate, 
Bonnell & Price (1995; sink particles), Watkins et al. (1995; 
Navier-Stokes viscosity), Morris & Monaghan (1997; time- 
varying artificial viscosity). Nelson & Langer (1997; scheme 
for treating the thermal physics), Owen et al. (1998; ten- 
sorial smoothing kernels), Inutsuka & Imaeda (2001; Go- 
dunov scheme for interparticle hydro forces) , Kessel-Deynet 
& Burkert (2000; scheme for following ionization fronts). 

The principle alternatives to SPH are finite difference 
(FD) codes, or hybrids which combine particles and cells. 
FD codes are much better at resolving large density and/or 
velocity gradients. They should probably be the method of 
choice for simulating flows with high Mach Number. How- 
ever, in three dimensions they are not as easy to implement 
as SPH, and they are not Lagrangian. Moreover, it should 
be born in mind that FD codes beneflt from having had very 
many more person-years of development than SPH. 

In FD codes, the high local resolution required for cal- 
culations of protostellar collapse can be obtained by means 
of nested grids. These nested grids can either be introduced 
at the outset to cover places where it is anticipated that ex- 
tra resolution will be required (e.g. Burkert & Bodenheimer 
1993); or they can be introduced and removed during the 
course of the simulation, whenever and wherever they are 
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required, as in Adaptive Mesh Refinement (AMR; Berger & 
Colella 1989, Truelove et al. 1997, 1998). With AMR, the 
resolution can, in principle, be increased indefinitely in re- 
gions where this is necessary. 

In standard SPH, the linear resolution is ~ {rn/p)^^^ , 
where m is the mass of an individual particle and p is the 
local density. In this paper we describe a rather straight- 
forward algorithm whereby the resolution of an SPH code 
can also, in principle, be increased indefinitely, by splitting 
particles. We test the method on the standard Boss & Bo- 
denheimer (1979) test, using the Jeans condition to trig- 
ger splitting. Particle splitting and merging were used by 
Meglicki, Wickramasinghe & Bicknell (1993) in their sim- 
ulations of accretion discs, and particle merging had been 
used earlier by Monaghan & Varnas (1988) in simulations 
of interstellar cloud complexes. However, in these earlier im- 
plementations splitting was used to maintain resolution in 
low-density regions, and merging was used to avoid follow- 
ing large numbers of particles in high-density regions (es- 
sentially the motivation for sink particles; Bate, Bonnell & 
Price, 1995). Therefore, the philosophy was very difi'erent 
from that adopted here, where high resolution is advocated 
in high-density regions to avoid violating the Jeans condition 
(i.e. to ensure resolution of the local Jeans mass). Moreover, 
no tests were reported in these earlier papers. 

Section 2 gives a brief review of standard SPH. Section 
3 sketches the implementation of Particle Splitting, at the 
microscopic level, and Section 4 describes how the best pa- 
rameters for Particle Splitting are determined. Sections 5 
and 6 describe two algorithms for Particle Splitting at the 
macroscopic level. In Nested Splitting (Section 5) all the 
particles in a prescribed sub-domain are split from some 
predetermined time igpUt onwards, so that in efi'ect a high- 
resolution simulation is performed inside the sub-domain us- 
ing initial and boundary conditions supplied by the coarser 
simulation in the overall computational domain. In On-The- 
Fly Splitting (Section 6), particles are split only when this 
is necessary according to some local criterion such as the 
Jeans Condition. On-The-Fly Splitting is the more compu- 
tationally efficient procedure, requiring no manual interven- 
tion and no prior knowledge of where high resolution will 
be required; it corresponds closely to AMR. Section 7 in- 
troduces the Jeans Condition, and Section 8 demonstrates 
how effectively and economically SPH performs the standard 
Boss & Bodenheimer (1979) rotating collapse test when Par- 
ticle Splitting is included. Section 9 discusses the results and 
compares them with those obtained using AMR. Our con- 
clusions are summarized in Section 10. 



2 BASIC SMOOTHED PARTICLE 

HYDRODYNAMICS AND TREE-CODE 
GRAVITY 

In SPH the evolution of the gas is predicted by following the 
motions of an ensemble of particles which act as dis- 

crete sampling points (e.g. Monaghan 1992). Each particle 
has associated with it a mass rui, smoothing length hi, posi- 
tion Vi, velocity Vi, and values for any other local intensive 
thermodynamic variables required to describe the state of 
the gas, for example density pi, pressure Pi, specific internal 
energy Ui, magnetic field Bi, etc. 



At an arbitrary position r, the value of any physical 
variable A can be obtained by interpolating over the values 
of A associated with nearby particles, using a smoothing 
kernel W , viz. 



The kernel is normalized, 

■oc 

W{s)A-KS^ds = 1 



(1) 



(2) 



and has compact support (specifically W{s > 2) = 0), so 
that the sum is over a finite number A/'ngib of neighbouring 
particles. 

The gradient of A can then be evaluated from 



V^(r) = ^ I 



Aiirii 
pihj 



W 



hi 



(3) 



where W'[s) = dW/ds. 

The motions of the particles are driven by interparti- 
cle forces representing the local pressure gradient, viscosity, 
gravity and magnetic field. Interparticle forces are formu- 
lated symmetrically so that linear and angular momentum 
are conserved. In our basic implementation we use the M4 
kernel introduced by Monaghan & Lattanzio (1985). The 
particles all have the same mass. We omit magnetic fields. 

The smoothing length hi is varied so that the kernel 
contains A/'jjgjb ~ 50 (±5) particles. This ensures that the 
spatial resolution becomes finer with increasing density. 

The density at the position of particle i is given by 



(4) 



where 



hij — 



hi + h-j 



ri 



The sum in Eqn. (j^) is over the ~ -^neih neighbours within 
2hij, and it includes particle i itself. 

In many applications, the pressure at particle i is given 
by a barotropic equation of state Pi = P{pi)- In this case 
there is no need to solve an energy equation. 

The gravitational acceleration of a particle can be ob- 
tained by a direct sum of the contributions from all the other 
particles. 



Grrij Ari, 



lAr, 



(5) 



However, there are two problems with this approach. The 
first problem arises because individual particles are very 
massive. As a result, close interactions between particles 
(small |Arij|) lead to artificially large accelerations which 
corrupt the simulation. This can be cured by softening the 
mutual gravity of any two particles having small separation. 
The form of softening we use assumes that the kernel de- 
scribes the distribution of a particle's mass, and then invokes 
Gauss's Gravitational Theorem. The effect is to introduce a 
term H^*(| Ar^j \/hij) into the sum of Eqn. (1) (see Eqn. (|), 
below), where 
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W*is) 



Ty(s')47rs'^ ds' 



(6) 



The equation of motion is integrated using a second- 
order Runge-Kutta scheme: 



The second problem arises because the cost of comput- 
ing the sum in Eqn. (^) for all the particles increases as 
*total (where i^otal total number of particles). Since 

the linear resolution, and the time-step, of an SPH simu- 
lation are proportional to i^^^^^^, the cost of increasing the 
resolution soon becomes prohibitive. An effective solution 
to this problem is to arrange the particles in an hierarchical 
octal tree structure, i.e. a nested grid of cells within cells. 
This is called Tree-Code Gravity (TCG; Barnes & Hut 1986, 
Hernquist 1987). 

The top cell of the tree is the entire computational do- 
main, assumed here to be a cube. The top cell is then parti- 
tioned into eight equal cubic sub-cells, these in turn are each 
divided into eight equal cubic sub-sub-cells, and so on. Parti- 
tioning is continued until the lowest cells in the tree contain 
either one particle or no particle. Then the mass Mi, centre 
of mass R/, linear size Li and quadrupole moment Q/ of 
each cell I are calculated and stored (along with some other 
cell parameters which are useful for finding neighbours). 

When the gravitational acceleration of particle i is be- 
ing computed, the contribution from a distant cell can be 
evaluated to sufficient accuracy using its M/,R/,Q/, and 
there is no need to consider its constituent particles individ- 
ually. For particle i at r^, cell / at position R/ is sufficiently 
distant to justify this approximation if 



Li < 



■^crit 



IRj 



(7) 



with S(,j.it ~ 0.576 (Salmon, Warren & Winckelmans 1994). 
The equation of motion for particle i is then 



IT 



E 



E 



2 ~^ 2 



|Ar,,P 



w 



hij 



Art 
|Ar, 

An 
lAr, 



(8) 



Here, the first sum represents short-range forces (hydrostatic 
and viscous forces) and is only over nearby particles. The 
second sum gives the gravitational acceleration, which is 
long-range, and is nominally over all particles. In practice, 
distant particles are grouped into clusters using the TCG 
algorithm, as described in the preceding paragraphs. 

Hij is the artificial viscosity term advocated by Lat- 
tanzio et al. (1985), viz. 



Hi 



where 



A-Vi 



^J.ij < 

l^ij > 



(Avij.Arij) hij 
|Ar,,|2 + 0.01^2. 



(9) 



(10) 



(ci -I- Cj)/2, and pij = {pi + pj)/2. We 



= a(rr,vr) ; 

n . n 

= ^^ + ai — ; 

( n+i n+i 



n + 1 



Vi -I- a- 



^Af : 



where r", v" and a" are the position, velocity and acceler- 
ation of particle i at time-step n. 

The code invokes a binary hierarchy of multiple particle 
time-steps, At\ = 2~*Afmax with k = 0, 1, 2, fcmax = 
30. For each particle i we evaluate a maximum time-step 
according to 



AU = 0.3 MIN. 



1 K 
IV.vL' WV 



hi 

NT 



1/2 



where 

cr, = c, + 1.2 {aci + /3MAX{/iij}) 



(11) 



(12) 



The particle is then given the largest time-step Atl from 
the binary hierarchy which is smaller than Ati, i.e. 



k 



Ati 



+ 1, 



AU 



Ati. 



(13) 



A more detailed description of the basic code can be 
found in Turner et al. (1995). 



3 PARTICLE SPLITTING 

The resolution of SPH is of order the local smoothing length 
h, and with A/'^gib ~ 50, h ~ (m/p)^^^. Therefore in order 
to improve the resolution of a simulation locally we need to 
reduce the masses of the individual particles, and hence to 
increase the overall number of particles. 

In Particle Splitting, this is done by replacing individual 
particles (parent particles) with small groups of particles 
(families of children particles). In the standard form of SPH, 
an individual particle is spherically symmetric, by virtue of 
having an isotropic smoothing kernel. Therefore we should 
distribute the family of children particles which replace it 
so that the sum of their smoothing kernels is approximately 
spherically symmetric. 

This is achieved by creating thirteen children, each with 



'"child = 



'"parent 
13 



(14) 



have used a = 1 and P — 1. 



The children are then placed on an hexagonal close-packed 
array, with the first child at the same position as the par- 
ent, and the other twelve children equidistant from the first; 
the method used to determine the optimum value for the 
distance £ between the first child and each of its siblings is 
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described in the next section. Finally the array of particles 
is rotated to an arbitrary orientation. 

The choice of 13 children is a compromise between two 
considerations. On the one hand, there should not be too 
large a disparity between parent mass and child mass, oth- 
erwise there will be greatly increased numerical diffusion 
wherever parents and children are neighbours. On the other 
hand, it is desirable that the collective density distribution 
of the family of children approximate to the spherically sym- 
metric density distribution of their parent. With 13 children 
there is a relatively small variance in the collective density 
of the children on a spherical surface centred on the position 
of the parent. 



4 FINE-TUNING THE FAMILY OF 
CHILDREN PARTICLES 

The children-particles are initially given smoothing lengths 

^child = (13)"^^^^parent ■ (15) 

so that they have the same net volume as their parent- 
particle. 

Two experiments have been preformed to determine the 
optimum value for I (the distance between the first child and 
each of its twelve siblings) . 

In the first experiment we take a microscopic ap- 
proach. We consider a single parent particle and compute 
the volume-weighted difference between its density and the 
combined densities of its children: 



all space 



Pparcnt 



(16) 



where the summation index i' represents the identifiers of 
the children. This quantity has its minimum value when I ~ 
1.9/i|,hjl(j, but it is still quite large, ~ 0.33mparent- 

In the second experiment we take a macroscopic ap- 
proach. 500 equal-mass particles are distributed randomly 
inside a cube. They are then settled using SPH with uni- 
form sound speed and periodic boundary conditions, until 
the density is very uniform, specifically (Jp/p < 0.01. Then 
we split all the particles simultaneously and re-settle them. 
For I — 1.5/i|,jj;[^, the value of Up/p is a minimum (~ 0.10) 
immediately after splitting, and the time required to settle 
back to Cp/p < 0.01 is also a minimum. The minimum of 
Cp/p immediately after splitting is quite shallow (it increases 
to ~ 0.12 for £ ~ 1.2/i(,jjjj^, and for £ ~ 2.0h^i^ii^), and this is 
largely testimony to the smoothing ability of the M4 kernel 
used with ^ 50 neighbours (Monaghan & Lattanzio 1985). 

We believe that this second experiment is the more rel- 
evant, since it is the macroscopic properties of the fluid 
which SPH seeks to mimic. In other words, it is the con- 
tinuous density and velocity fields (which at most positions 
are compounded by contributions from ~ A^neib particles) 
that SPH must aspire to reproduce, and not the density of 
a single parent-particle. 

All subsequent tests are therefore performed with £ = 
1-5'ichild- 



5 NESTED SPLITTING 

One way to implement Particle Splitting is to start a stan- 
dard simulation at f = 0, then a,t t — tsplit(— 0) ^'^ identify a 
part of the overall computational domain, i.e. a sub-domain, 
where improved resolution is required or desired, and to split 
all particles in this sub-domain. In the sequel we refer to the 
initial standard simulation and the particles in it as - re- 
spectively - the coarse simulation and the coarse particles. 
We refer to the simulation after fgpiit and the split particles 
in the sub-domain as - respectively - the fine simulation and 
the fine particles. 

At time t = tspm the newly-created fine child-particle 
i' is given a velocity by summing contributions from the 
~ 50 neighbours (index j) of - and including - the coarse 
parent-particle (index i): 



W 



(17) 



Thereafter {t > ispiit)i any coarse parent-particle which en- 
ters the sub-domain is immediately split and the resulting 
fine child-particles are again given velocities by summing 
over the neighbours of the coarse parent-particle. The den- 
sities and accelerations of fine particles are calculated using 
the standard SPH equations (Eqns. (|) & (|)). 

Near the boundary of the sub-domain, there are coarse 
particles which have fine particles as neighbours, and vice 
versa. Because of this, we have to revise the way in which 
the smoothing lengths of particles are adjusted. The results 
are greatly improved if, instead of requiring the number of 
neighbours of particle i to be ~ 50, we require the total mass 
of the neighbours of particle i to be ~ 50 times the mass of 
particle i, i.e. 



{mj} ~ 50mi 



(18) 



For fine particles well inside the sub-domain, and for coarse 
particles well outside the sub-domain, this makes no differ- 
ence, but it is important for particles near the boundary, 
including coarse particles which are about to be split. 



6 ON-THE-FLY SPLITTING 

A second way to implement Particle Splitting is to iden- 
tify particles whose resolution, ~ {m/pY^^, is insufficient 
for the problem in hand, and split them on-the-fiy. To do 
this, we must recognize what the deterministic physics of 
the problem is, and quantify - as a local function of state 
- the minimum resolution required to capture this physics. 
Then, whenever the resolution of a particle approaches the 
local minimum resolution, that particle is split. 

As in Nested Splitting, the velocities of the new 
children-particles are evaluated by summing contributions 
from the neighbours of the parent particle. Successive gener- 
ations of splitting are possible, yielding particles with masses 
m/13, m/13^, m/13'', etc. and hence linear resolution im- 
proved by factors 2.4, 5.5, 13, etc. 

As in Nested Splitting, better smoothing is obtained if 
the desired number of neighbours is evaluated according to 
Eqn. (§1). 

In the context of star formation, gravitational collapse 
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and fragmentation are crucial physical processes, and there- 
fore the code must always be able to resolve the Jeans mass. 
This is called the Jeans Condition (Truelove et al. 1997, Bate 
& Burkert 1977). Since the Jeans mass is a local function 
of state, it is straightforward to formulate the Jeans Condi- 
tion, and to split particles when they are about to violate it. 
We explore On-The-Fly Particle Splitting triggered by the 
Jeans Condition, in the next two sections. 



7 THE JEANS CONDITION 

Truelove et al. (1997) have shown that simulations of 
protostellar collapse and fragmentation performed using 
AMR only converge when the local Jeans length Ajgg^^g ~ 
a{Gp)~^^^ is resolved (where a is the local sound speed). 
Specifically the linear resolution Ar must everywhere sat- 
isfy 



(19) 

Otherwise artificial fragmentation can occur and/or real 
fragmentation can be suppressed. Eqn. ( p^ is the formu- 
lation of the Jeans Condition appropriate for FD codes. 

Bate & Burkert (1997) have shown that there is a sim- 
ilar Jeans Condition for SPH simulations: the minimum re- 
solvable mass ~ 2A/'neibm must everywhere be less than the 
local Jeans mass 



G-3/V-^/'a^ or 



G^PiVat 



< 



1 



(2a; 



10" 



neib/ 



(20) 



Otherwise artificial fragmentation may occur, and real frag- 
mentation will definitely be suppressed. 

Whitworth (1998) has shown that provided the smooth- 
ing kernel is sufficiently centrally peaked, SPH simulations 
of collapse are unlikely to suffer from artificial fragmenta- 
tion. However, it remains the case that Eqn. (||) must be 
satisfied if real fragmentation is to be modelled. 



8 ROTATING COLLAPSE TEST 

A standard test of star formation codes is the one proposed 
by Boss & Bodenheimer (1979; hereafter BB79). In this test 
the initial conditions are a uniform-density isothermal cloud 
having mass A/0, radius 0.02 pc, density 2 x 10~^* g cm""^ 
and isothermal sound speed 0.17 km s^^ (so the ratio of 
thermal to gravitational energy is ~0.26). The cloud then 
has an m = 2 azimuthal density perturbation with 10% am- 
plitude imposed on it, and it is set to rotate at uniform 
angular speed 7.2 x 10~^^ rad s~^ (so the ratio of rotational 
to gravitational energy is ~0.16). There is no external pres- 
sure, and the gas remains isothermal during the subsequent 
evolution, until the density rises above a critical value p^rit' 
when it switches to being adiabatic and therefore heats up. 

Klein et al. (1999) have performed this test, with Pcrit = 
5 X 10"^'* g cm~^, using AMR, and shown that the cloud 
should collapse to produce an elongated structure. This 
structure then evolves into a binary system with a bar be- 
tween the components. As long as the gas remains isother- 
mal, the binary components and the inter-connecting bar 
condense into filamentary singularities, and the bar does not 
fragment (Truelove et al. 1998) . These results agree with the 



theoretical analysis of Inutsuka & Miyama (1992), but con- 
trast with earlier simulations of the BB79 test using both 
FD and SPH codes. In particular, the simulations reported 
by Burkert & Bodenheimer (1993) resulted in the bar frag- 
menting. The reason for this appears to have been that their 
simulations did not satisfy the Jeans Condition. 

Bate & Burkert (1997) have subsequently repeated the 
BB79 test with p^rit = lO"^"* g cm"^, using SPH and 80,000 
particles, sufficient to satisfy the Jeans Condition up to a 
density ~ Pcrit- They now find essentially the same result 
as Truelove et al. (1998) and Klein et al. (1999), confirming 
that the earlier results of Burkert & Bodenheimer (1993) 
were compromised by violating the Jeans condition. Thus 
the BB79 test appears to be a stringent test of whether a 
code is satisfying the Jeans Condition. 

We have performed three SPH simulations of the BB79 
test. We have required the gas to remain isothermal up to a 
density p^rit = 5 X 10"^'^ g cm~^. This is two orders of mag- 
nitude higher than the value used by Klein et al. (1999), and 
so in our simulations the Jeans mass and length fall to val- 
ues 10 times smaller, i.e. Afjgg^^g ~ 1O~*M0 and i? jeans ^ 2 
AU. Thus we are subjecting our code to an even sterner test 
than that advocated by Klein et al. (1999). If one accepts 
the notion of opacity-limited fragmentation, then in nature 
the Jeans mass should not fall below ~ IO^^A/q. Therefore 
- in this limited sense - we are also subjecting our code to 
a sterner test than that set by nature. 

In the standard simulation there are 600,000 particles 
from the outset, sufficient to satisfy the Jeans condition 
without Particle Splitting. In the other two simulations there 
are only 45,000 particles at the outset, and the Jeans con- 
dition is accommodated by implementing particle splitting, 
first Nested Particle Splitting and then On-The-Fly Parti- 
cle Splitting. The Nested Particle Splitting simulation ends 
up with ~ 320,000 particles, and the On-The-Fly Particle 
Splitting simulation ends up with ~ 140,000 particles. 

We follow the suggestion of Truelove et al. (1997, 1998) 
that numerical results should be compared, not at precisely 
the same elapsed time, but instead when the peak density in 
the computational domain, Ppeaki reaches the same value. In 
particular we present results when Ppeak — Pcrit = 5 x 10^^^ 
g cm~^, and when Ppeak — 2 X 10~^ g cm""^. The reason for 
this is that the evolution of a simulation appears to progress 
slightly more rapidly if the resolution is higher. The same 
effect was reported by Truelove et al. (1997, 1998) for their 
AMR simulations. 

The plots are all grey-scale column-density images 
through the central 0.004 pc x 0.004 pc (800 AU by 800 
AU) of the computational domain, viewed down the rotation 
axis. The calibration of the grey-scale is given in the figure 
captions. Times are given in terms of the initial freefall time 
tpF ^ 0.03 Myr. 



8.1 Standard high-resolution SPH simulation 
(no Particle Splitting) 



In Figure 1 we show the results of a simulation of the BB79 
test with pcj-it = 5 x 10~^^ g cm~^, performed using our 
standard SPH code with 600,000 constant-mass particles. 
This simulation has sufficient particles to satisfy the Jeans 
condition throughout the simulation without Particle Split- 
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Figure 1. The BB79 test performed using standard SPH with 600,000 constant-mass particles and anticlockwise rotation: column- 
density images of the central 0.004 X 0.004 pc^ of the computational domain, looking down the rotation axis, (a) At time t = 1.237tFF, 
when Ppga^jj — Pcrit = 5 X 10"^'^ g cm~^ ; the grey-scale is logarithmic, with sixteen equal intervals from 3.60 X 10^'^ H2 cm^^ to 
4.56 X 10^5 H2 cm ^ (or equivalently, adopting solar composition, 1.44 gem ^ to 1.82 X lO'^gcm (b) At time t = 1.245tFF when 
/'peak — 2 X lO^'^* g cm~'^ ; the grey-scale is logarithmic, with sixteen equal intervals from 3.72 X lO'^^ H2 cm"'^ to 1.55 X lO^'^ H2 cm~^ 
(or equivalently 1.49gcm~^ to 6.19 X lO^gcm"^). 



ting. We note that our SPH code has been developed entirely 
independently of that used by Bate & Burkert (1997) and 
differs from that code in several major regards, in particular 
the integration scheme and the gravity solver. 

Figure la shows a grey-scale column- density image of 
the centre of the computational domain at the moment the 
density in the binary components Ppeak passes pcriti i-^- 
a,t t = 1.237tFF. Both the binary components and the bar 
connecting them approximate to filamentary singularities, 
and there is no sign of bar fragmentation. 

Figure lb shows the same region at the end of the sim- 
ulation, t = 1.245tFF. By this stage the gas in the densest 
parts, i.e. the binary components, has become adiabatic and 
heated up. As a result the binary components are approxi- 
mately circular in projection. However, the interconnecting 
bar is still isothermal; it continues to approximate to a fila- 
mentary singularity and shows no sign of fragmenting. 

We re-iterate that this simulation extends the isother- 
mal evolution to p^rit = 5 X 10~^^ g cm~^, which is two 
orders of magnitude higher than the critical density used by 
Klein et al. (1999). Therefore it constitutes a very stern test 
of the code's validity. 

8.2 SPH simulation with Nested Particle Splitting 

Figure 2 shows the result of a simulation of the BB79 test, 
again with pcrit = 5 x 10"^'^ g cm~'^, but now performed 
with only 45,000 equal-mass particles at the outset. Nested 
Particle Splitting is applied after tspia = 1.244tFF, within 
a cylindrical sub-domain of radius 0.003 pc (600 AU) and 



height 0.003 pc (600 AU). The symmetry axis of the cylinder 
is aligned with the rotation axis of the cloud. Its size is 
chosen so as to contain all regions which become so dense 
that without Particle Splitting they would eventually violate 
the Jeans Condition, i.e. the final binary components, which 
have an orbit of radius 0.002 pc (400 AU). 

Figures 2a and 2b correspond to times t — 1.258tFF 
(when Ppeak — Pcrit) ^^'^ ^ ~ 1.265fpF (the end of the sim- 
ulation at Ppeak — 2 X 10~^ g cm~'^). By comparing Figs. 1 
and 2, we see that the simulation with Nested Particle Split- 
ting reproduces the crucial features of the evolution found in 
the standard high-resolution simulation described in Section 
8.1. As before, a binary system forms with a bar between the 
components, and as long as the gas remains isothermal the 
binary components and the bar evolve towards filamentary 
singularities, with no tendency for additional fragments to 
condense out of the bar. There is a small density peak at 
the centre of the bar in the first frame (2a; Ppg^k — Pcrit )> 
but this quickly disperses, and therefore we do not consider 
it to be critical. 

By the end, approximately half the particles have been 
split, so there are ~ 320,000 particles in total. The simu- 
lation with Nested Particle Splitting therefore requires less 
than half the memory and approximately 25% of the CPU 
used by the standard simulation (with 600,000 particles, 
§8.1). 
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Figure 2. As Fig. 1, but performed using only 45,000 particles initially and then Nested Particle Splitting in a central cylindrical 
sub-domain of radius 0.003 pc and height 0.003 pc. (a) At time t = 1.258tFF, when Pp^gjj^ — Pcrit = 5 X 10~^^ g cm"'' . (b) At time 
t = 1.265tFF when Ppgak - 2 X 10"^ g cm-^ . 




Figure 3. As Fig. 1, but performed using only 45,000 particles initially and then On-The-Fly Splitting triggered above p = p^piit = 
3 X lO"!"" g cm~^. (a) At time t = 1.259iFF, when Ppg^k — Pcrit = 5 X 10~^^ g cm"^ . (b) At time t = 1.270tFF when Ppgak — 3 x 10~^ 
g cm~^ . 
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8.3 SPH simulation with On-The-FIy Particle 
Splitting 

Finally Fig. 3 shows the result of a simulation of the BB79 
test, again with p^rit = 5 X 10"^'^ g cm~'^, and again per- 
formed using only 45,000 equal-mass particles at the outset, 
but now with Particle Splitting applied On-The-Fly, in re- 
sponse to imminent violation of the Jeans Condition. 

Figs. 3a and 3b correspond to times t — 1.259tFF (when 
Ppeak — Pcrit) ^''id t = 1.270tFF (the end of the simulation 
at Ppeak = 3 X 10~^ g cm~^; a slightly higher final density 
than in the other two simulations, see below). Again the 
simulation reproduces all the critical features reported by 
Truelove et al. (1998) and Klein et al. (1999), in particular 
no tendency for additional fragments to condense out of the 
bar. Once more there is a small density peak at the centre 
of the bar in the first frame (3a; Ppg^k — Pcrit)' which then 
quickly disperses. 

The only significant difference is a fanning-out of the 
bar close to the binary components, right at the end of the 
simulation. This is a real effect, due to the tidal shearing of 
the bar by the inspiralling binary components. It only ap- 
pears in this simulation because the evolution has been fol- 
lowed to somewhat higher density (3 x 10~®gcm~'^ instead 
of 2 X 10~^ gcm~^); hence the binary components have got- 
ten further out of alignment with the bar, where they can 
exert a disruptive shear on the ends of the bar. 

At the end, less than one fifth of the original particles 
have been split, so there are ~ 140,000 particles in total. 
The simulation with On-The-Fly Splitting therefore requires 
less than one quarter the memory and less than 10% of the 
CPU used by the standard simulation (with 600,000 parti- 
cles, §8.1). 



9 DISCUSSION 

It appears that Particle Splitting is a viable option for in- 
creasing the resolution, locally, during an SPH simulation of 
self-gravitating collapse. In this context, the most stringent 
and apposite test of the algorithm is the BB79 test, where 
extra resolution is required to avoid violating the Jeans Con- 
dition (and hence to avoid artificial fragmentation). Accept- 
able results are obtained for the BB79 test, even when the 
gas is programmed to stay isothermal up to Pcrit = 5 x lO"^'^ 
g cm~^. Perturbations due to the splitting process are tran- 
sient, and are effectively and quickly damped. Moreover, the 
simulation is accomplished with much less memory and CPU 
than a standard simulation, particularly with the On-The- 
Fly implementation. We note also that the savings in mem- 
ory and CPU are likely to be far greater in realistic simula- 
tions of star formation, where only a fraction of the matter 
in a cloud collapses to form stars, and therefore a far smaller 
fraction of particles needs to be split and to be followed with 
a small time-step. 

Other tests have been performed using the two Parti- 
cle Splitting algorithms. In particular we note the following, 
(i) We have evolved a static, stable isothermal sphere, con- 
tained by an external pressure Pexti ^-nd shown (a) that it 
remains stable, and (b) that its density profile is not cor- 
rupted. The configuration treated has mass Mq, radius 0.01 
pc, and temperature 7.9 K. It is quite centrally condensed. 



with Pcentral - 3Pedge and Pext ==; 3 X 10 erg cm ^ (ii) 
We have also followed the collapse of an isothermal cloud 
which initially has uniform density, and quickly evolves to- 
wards an inverse-square density profile (p (x r~^). The cloud 
has mass Mq and temperature 7.9 K. Initially it has radius 
0.016 pc and uniform density ~ 4 x 10~^^ g cm~^. The effect 
of Particle Splitting on the ensuing collapse is negligible. 

Particle Splitting might be applied in other situations 
where increased resolution is required locally. The Nested 
implementation requires prior knowledge of where extra res- 
olution will be required. The On-The-Fly implementation re- 
quires that the condition for requiring extra resolution can 
be formulated as a local function of state. However, more 
development and testing is needed to establish the stability 
of the method in these more general situations. 

One difference between Particle Splitting and AMR is 
that Particle Splitting - at least in its present form - is irre- 
versible. If fiuid fiows into, and then out of, a region where 
extra resolution is required, particles get split and stay split. 
However, there is no a priori reason why particles should 
not be merged in regions where they are delivering much 
higher resolution than is required by the local physics, and 
indeed this was done successfully by both Monaghan & Var- 
nas (1988) and by Meglicki, Wickramasinghe and Bicknell 
(1993), using a grid. We have not yet attempted to do this, 
so we cannot comment further on its feasibility. 

The concern is that, without merging, if a single low- 
mass particle is surrounded by high-mass particles, Eqn. ( p^ 
dictates that it has very few neighbours (typically only 4 if 
the surrounding particles are 13 times more massive) . There- 
fore the low-mass particle might have rather noisy density, 
pressure, etc. because of small-number statistics and sud- 
den changes in its neighbour list. However, this will be to a 
large extent offset by two factors. First, by implication the 
fiuid is well resolved by the large particles, so they have a 
well relaxed distribution and the neighbourhood of the low- 
mass particle should be quite constant. Second, as regards 
the macroscopic continuum properties of the fiuid, within 
the sphere of infiuence of the low-mass particle (i.e. within 
its smoothing kernel) it makes only a ~ 2% contribution to 
the local thermodynamic variables, and so fiuctuations in its 
individual properties are very diluted. 

Shocks and ionization fronts can be treated effectively 
using ellipsoidal kernels (Owen et al. 1998; Francis 2001). 
Moreover, in the context of star formation, irreversibility 
may not be a serious problem, because once gravitational 
collapse and fragmentation are underway the local Jeans 
mass is unlikely to increase much, and so the resolution re- 
quired by the Jeans Condition is unlikely to decrease. We 
are currently developing a version of our code which treats 
magnetic fields and ambipolar diffusion (Hosking & Whit- 
worth, in preparation), but until this has been fully tested 
without particle splitting it would not be sensible to combine 
the two. 



10 CONCLUSIONS 

We have described how Particle Splitting can be imple- 
mented in SPH, and we have shown that it offers a stable 
and economic way to increase the local resolution of an SPH 
simulation of self-gravitating collapse, and hence avoid vio- 
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lating the Jeans Condition. This makes Particle Splitting a 
useful refinement for SPH simulations of gravitational col- 
lapse and fragmentation. 
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